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Abstract 

a 

A well established framework of an uncoupled hierarchical modeling approach is 
adopted here for the prediction of macroscopic material parameters of the General- 
ized Leonov (GL) constitutive model intended for the analysis of flexible pavements 
at both moderate and elevated temperature regimes. To that end, a recently intro- 
duced concept of a statistically equivalent periodic unit cell (SEPUC) is addressed 
to reflect a real microstructure of Mastic Asphalt mixtures (MAm). While mastic 
properties are derived from an extensive experimental program, the macroscopic 
properties of MAm are fitted to virtual numerical experiments performed on the 
basis of first order homogenization scheme. To enhance feasibility of the solution of 
the underlying nonlinear problem a two-step homogenization procedure is proposed. 
Here, the effective material properties are first found for a mortar phase, a compos- 
ite consisting of a mastic matrix and a fraction of small aggregates. These properties 
are then introduced in place of the matrix in actual unit cells to give estimates of 
the model parameters on macroscale. Comparison with the Mori-Tanaka predictions 
is also provided suggesting limitations of classical micromechanical models. 
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1 Introduction 



As seen in Fig. [T](a), asphalt mixtures represent in general a highly heteroge- 
neous material with complex microstructure consisting at minimum of mastic 
binder, aggregates and voids. When limiting our attention to Mastic Asphalt 
mixtures, used typically in traffic arteries of substantial importance, the frac- 
tion of voids becomes negligible. A binary image of such a two-phase material 
system plotted in Fig. [TJb) is then readily available. The literature offers sev- 
eral distinct routes taking advantage of such a representation. 




(a) (b) (c) 

Fig. 1. (a) A real microstructure of an asphalt mixture, (b) Original binary image, 
(c) Improved binary image 

Phenomenological constitutive models enhanced by introducing a microstruc- 
ture linked aggregate distribution tensor [Tf2] were proposed, following the 
footsteps of isotropic viscoelastic or viscoplastic models developed since the 
early 1990s [3 4 5,6,7, to cite a few], to account for an intrinsic aggregate struc- 
ture anisotropy. Recently, detailed micromechanical models developed on the 
basis of image processing of either two-dimensional digital photography or 
three-dimensional X-ray computed tomography have been favored for the per- 
formance assessment of asphalt concrete composites. Typically, a sufficiently 
large sample of microstructure of a two-phase composite system, transformed 
into a binary image, was adopted to serve as a representative volume ele- 
ment (RVE). Either direct application of these models [8]f9]fTD] or application 
of simplified artificial microstructures of the same type [TTJ was examined to 
provide estimates of the macroscopic response. While capturing both the local 
and overall behavior sufficiently accurately, most of these models suffer from 
computational inefficiently. Exploitation of the true potential of classical aver- 
aging micromechanical schemes combined with "bottom-up" hierarchical ho- 
mogenization technique thus appears rather natural. As an example we recall 
a successful application of the Mori-Tanaka method in upscaling viscoelastic 
properties of asphalt mixtures at low temperatures [T2ll3M14|[T5] . 

Growing from the roots planted by the authors in the early 2000s [T6lll7|ll8j the 
present contribution combines these various aspects of micromechanical mod- 
eling into a relatively simple, yet reliable and efficient computational scheme. 
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While still incorporating the prominent "bottom-up" uncoupled multiscale ho- 
mogenization scheme, the effective properties on individual scales are found as 
volume averages of local stress and strain fields developed inside the so-called 
Statistically Equivalent Periodic Unit Cell [T9] . 



Rendering the desired macroscopic constitutive model that describes the ho- 
mogenized response of M Am to general loading actions thus endeavors to the 
formulation of a suitable micromechanical model on individual scales and to 
associated experimental work being jointly the building blocks of the upscal- 
ing procedure. Three particular scales shown in Fig. [2] are considered in the 
present study: 




Fig. 2. Three distinct scales of Mastic Asphalt mixture 

Although the mastic-phase itself is a composite consisting of a filler and 
a bituminous binder, it is assumed in the present study to be well repre- 
sented by a temperature and rate dependent homogeneous isotropic mate- 
rial. Since limiting our attention to elevated temperatures exceeding 30°C, 



the GL model discussed in Section 3.1 is exploited to provide for experi- 



mentally observed nonlinear viscoelastic behavior of bituminous matrices. 
The required experimental program to calibrate the model parameters is 
outlined in Section 13.21 



As presented in Section 2.1 , the mortar-phase naturally arises through the 
process of removing aggregates from original microstructure, Fig.[l](b), pass- 
ing 2.26 mm sieve. Introducing the principal assumption of the adopted 
upscaling procedure - the GL model is suitable for governing the homog- 
enized response on every scale - allows us to derive the model parameters 
on the mortar-scale by running a set of virtual numerical experiments on 
the simplified periodic hexagonal array model plotted in Fig. [2j The selec- 
tion of this geometrical representation of the mortar composite is purely 
an assumption building upon the conclusion that at least for the mastic- 
phase and low-temperature creep the response is independent of the filler 
shape and mineralogy [13J. Various computational strategies delivering the 



searched macroscopic response are presented in Sections 2.2.1 - 2.2.3 
The SEPUC, also seen in Fig. [2j is selected to predict the macroscopic re- 
sponse of MAm, where large aggregates are bonded to a mortar-phase being 
homogeneous and isotropic. The elliptical shape of aggregates has already 
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been successfully used in detailed micromechanical simulations presented 
in [TT]. The issue of constructing SEPUC for the class of Mastic Asphalt 



mixtures is briefly addressed in Section 2.3 



The proposed upscaling procedure consistent with Fig. is then applied in 



Sections |3.3| and |3.4| to deliver the necessary parameters of the homogenized 
GL model. 



In the following text, a, a and A denote a vector, a symmetric second-order 
and a fourth-order tensor, respectively. The standard summation notation is 
adopted, i.e., A : B and A : b denote the sums A^B^s and A^^m, a : b 
being aijbij represents a scalar quantity and a ■ b stands for ciijbj, where the 
summation with respect to repeated indices is used. The symbol {a} is reserved 
for a column matrix or a vectorial representation of symmetric second-order 
tensor while the notation [L] is employed for a matrix representation of a 
fourth-order tensor [20J. 



2 Concept of statistically equivalent periodic unit cell 



It has been demonstrated in our previous work, see e.g. [T7]l21f22j that image 
analysis of real, rather then artificial, material systems plays an essential role 
in the derivation of a reliable and accurate computational model. This issue is 
revisited here to get, although considerably simplified, a reasonably accurate 
computational model (RVE^SEPUC) representing the real microstructure of 
Mastic Asphalt mixtures. 



2.1 From real microstructure to computationally feasible RVE 



Notice a 1,600x1,600 pixel resolution binary image in Fig.[T](b) created through 
the process of color segmentation of the original image in Fig. [T|(b) using the 
freeware graphical software GIMP. Visible flaws, e.g. defects inside large ag- 
gregates, required, however, further modification. To that end, an automated 
numerical tool was developed to determine a boundary of each aggregate, fill 
defects inside large aggregates and separate all aggregates as much as possible. 
The resulting improved image seen in Fig. [l](c) then allows for microstructure 
quantification in terms of various statistical descriptors. 

To open this subject, we consider first cumulative distribution functions of 
aggregates plotted in Fig. |3j The two graphs suggest a relatively large amount 
of small stones (cumulative distribution function of number of aggregates) 
which might seem, at least from their volume fraction point of view (cumula- 
tive distribution function of volume fraction of aggregates), almost negligible. 
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Fig. 3. Cumulative distribution function 

Eliminating the small aggregates to simplify the original microstructure is nat- 
urally promoted. Several examples of such microstructures appear in Fig. |4| 
If admitting removal of all fragments smaller than 1,200 pixels in area (ag- 
gregates passing 2.36mm sieve) we reduce the total number of stones by 97% 
while the original volume fraction of aggregates dropped down only by 21%. 
However, neglecting the small aggregates completely would severely underesti- 
mate the final macroscopic predictions. A two-scale homogenization approach 
is therefore inevitable. 




(c) (d) 

Fig. 4. Examples of binary images of original microstructure after eliminating stone 
fragments smaller than (in area): (a) 150px, (b) 300px, (c) 600px, (d) 1200px 
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2.2 Upscaling material properties towards Mastic Asphalt Mixes 



Suppose that the micro-structure in Fig. W[d) is sufficient to provide reason- 
ably accurate predictions of macroscopic response of MAm. With reference to 
Fig. [2] and to findings presented in [13] we further assume the mortar phase 
to be well represented by a periodic hexagonal arrangement of aggregates of a 
circular cross-section. The resulting Mastic Asphalt mixture then consists of 
41% volume of large aggregates and 59% volume of mortar phase, whereas the 
mortar-mastic composite was calculated to have 34% volume of small stones 
and 64% volume of mastic phase. The latter composite is shown in Fig. |5]^d). 
For further reference the subscript s will be reserved for aggregates (stones) 
whereas the subscript m will be used to denote the binder (either mortar or 
mastic matrix). 





(a) A s = 0.19A (b) A s = 0.25A (c) A s = 0.31A (d) A s = 0.34A 

Fig. 5. PUC of mortar phase containing stone fragments smaller than (in area): (a) 
150px, (b) 300px, (c) 600px, (d) 1200px 

Given the geometrical models for individual scales now opens the way to the 
derivation of the effective elastic properties of MAm employing the uncou- 
pled "bottom-up" upscaling scheme. While confirming applicability of this 
approach, these results will further serve as a point of departure for generat- 
ing the statistically equivalent periodic unit cells and at the same time will 
provide a source of data to check the quality of the results obtained for these 



artificial microstructures, see Section 2.3 



A number of homogenization techniques is available in the literature. Three 
particular representatives will be now briefly reviewed. Owing to the solu- 
tion of the underlying nonlinear viscoelastic problem, the formulation will be 
presented in an incremental format. 



2.2.1 The Mori-Tanaka method 

Consider a two-phase composite medium Q loaded on its outer boundary S 
by an affine displacement field Auq = AE ■ x consistent with the macroscopic 
uniform strain increment AE. The local constitutive equation can be written 

as 

A<r(x) = L(x) : Ae(x) + AA(x) or Acr r = L r : Ae r + AA r , (1) 
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when limiting our attention to a piecewise uniform variation of the stress and 
strain fields within individual phases; L r is then the stiffness tensor of the phase 
r = s,m and the increment of eigenstress AA r represents contributions caused 
by various non-mechanical physical sources. In the present study we admit only 
the creep strains developed in the binder phase thus leaving the aggregates 
to remain linearly elastic. In the framework of Dvorak's Transformation Field 
Analysis [23.24] the local strain increments then read 



Ae s = A S : AE + D sm : Afi m , (2) 
As rn A m . AE + D mm : A^ m , (3) 

where A r ,D rm ,r = s,m, are the mechanical strain localization tensors and 
transformation influence functions, respectively. The phase eigenstrain fi m is 
introduced to account for a nonlinear viscoelastic behavior of the binder. Note 
that for a two-phase medium, the transformation influence functions are di- 
rectly available in terms of the strain localization tensors [23] . 



Next, writing the macroscopic constitutive law as 

AS = L hom : AE + AA, (4) 

and realizing that 

AS = c s Acr s + c m Acr m , (5) 
yields the macroscopic stiffness matrix and the macroscopic eigenstress [23] as 

L c s L s . A s -|- c m L m . A m , (6) 
AA [c s (L s . D sm . M m ) -|- c m (L m . D mm . M m )] . AA m -I- c m AA m , 

v. / 

=0 

Cm,A m C m \- m /J, Jn . (7) 

Herein, the strain localization tensors A r are found from the Mori-Tanaka 
(MT) method [25]. In Benveniste's reformulation [26], the method approxi- 
mates the particle interactions by loading an isolated heterogeneity bonded 
to an unbounded matrix by a uniform stress or strain found in the matrix. 
Introducing the partial strain concentration (localization) tensor T s for the 
phase s (note that T m = I, the identity tensor), the local strain increment in 
aggregates is, in the absence of inelastic effects, provided by 

Ae s = T s : Ae m . (8) 

Next, writing AE = c s Ae s + c m As m readily provides 

As m = [c m l + c s T s ]" 1 : AE, (9) 
Ae s = (T s : A m ) : AE, (10) 
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which gives the strain localization tensors in the form 

A m [c m l ~\~ c s T s ] A s T s . A m . (1-0 

Closed form solutions for the localization tensor T s are provided in the litera- 
ture for certain types of heterogeneities, see e.g. [22] • Here, the model of aligned 
infinite cylinders of circular cross-sections, consistent with the assumed plane 
strain conditions, is adopted. Since used later in the Appendix, we present 
here the homogenized in-plane shear modulus m in the form 

^ r m>s m m(k m + 2m m ) + k m m m {c s Tfi s + c m m m ) (12) 

kmX^rn \k"m d" 2w m ) (c s 77l m -|- C m T7i s ) 

where m r ,r = s,m are corresponding shear moduli of individual phases. If 
both phases are also isotropic we get k r = m r /(l — 2z/ r ), where v r is the Poisson 
ratio. Finally note that all localization and transformation tensors mentioned 
in this section are functions of the instantaneous material properties of the 
binder phase, the current value of the viscoelastic shear modulus in particular. 



2.2.2 The Fast Fourier Transform method 

The Fast Fourier Transform (FFT) method proposed in [28] can be imagined 
as a passage between classical micromechanics models and the concept of 
periodic homogenization. To introduce this approach consider again a two- 
phase composite medium subjected to the same boundary conditions as in the 
previous section. 




General problem Step I Step II 



Fig. 6. Body with prescribed surface displacements including eigenstresses 

As suggested by Hashin and Shtrikman the local stress and strain fields 
in Eq. can be found from two auxiliary boundary value problems, Fig. |6} 
The procedure starts by assuming a geometrically identical body with a certain 
reference homogeneous, but generally anisotropic, medium l_o and the same 
prescribed displacements. The corresponding uniform strain E and stress S 
fields are related through constitutive law as 

AS = L : AE in ft, Au = Au on S. (13) 
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Following the Hashin-Shtrikman idea, we introduce the symmetric stress po- 
larization tensor r such that 



Acr(x) = L : Ae(x) + Ar(x). 

In addition, denote 

u* = u — uq in Q, u* = on S, e* = e — E, in Q. 
The increment of stress polarization tensor r then becomes 
At(x) = (L(x) - L ) : Ae(x) - AA(x). 



(14) 



(15) 



(16) 



If the polarization stress is known, the fluctuation part of the strain field e* 
can be obtained via Green's function l~° for a given reference medium in the 
form, see e.g. 



Ae*(x) = - / r°(x - x') : Ar(x') dx'. 
Jn 



(17) 



Combining Eqs. (15)3,(16) and (17) we obtain the so called periodic Lippmann 
Schwinger integral equation for a given reference medium as 



Ae(x) + / 



r°(x 



[(L(x') 



Ae(x') - AA(x)] dx' = AE. (tt 



This equation can be solved by the following iterative procedure 



As 



k+l, 



AE 



r°(x 



'L(x') - L°) : Ae k W 



A\ k W) 



dx'. 
(19) 

Typically, the Fast Fourier Transform (or rather its discrete version when 
dealing directly with binary images) is employed to solve the above equa- 
tion. Details on actual numerical implementation can be found in (28]. In the 
absence of inelastic effects the above equation simplifies as 



£ fe+1 (x) 



E 



r°(x-x'): [(L(x') 



s fc (x': 



dx'. 



(20) 



For plane-strain problem the homogenized elastic stiffness matrix |_ hom can be 
found from the solution of four successive elasticity problems. To that end, 
the RVE is loaded, in turn, by each of the four components of E, while the 
other three components vanish. The volume stress averages normalized with 
respect to E then furnish individual columns of |_ hom . 



2.2.3 First order homogenization of periodic fields 

Consider a material representative volume element now defined in terms of a 
periodic unit cell (PUC). Similar to the previous sections, suppose that the 
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PUC is subjected to boundary displacements Awo resulting in a uniform strain 
AE throughout the body. The strain and displacement fields in the PUC then 
admit the following decomposition, recall also Eq. (15)3, 



Ati(x) = AE • x + Au*(x), 
Ae(x) = AE + Ae*(x), 



Vxefi, m = m Vx6S, 

Vx g n. 



(21) 
(22) 



The first term in Eq. (21) corresponds to a displacement field in an effective 



homogeneous medium with the same overall properties as the composite. The 



fluctuation part u* enters Eq. (21) as a consequence of the presence of het- 

for further 



erogeneities and has to disappear upon volume averaging, see 
discussion. This condition is met for any periodic displacement field with the 
period equal to the size of the unit cell under consideration, (32J and references 
therein]. The periodicity of u* further implies that the average of e* in the 
unit cell vanishes as well. Hence 



(e(x)> = £J+( e *(x)> 



(£*(*)> 



ttJn 



e*(x) dx = 0. 



(23) 



Derivation of the macroscopic response then relies on the application of Hill's 
lemma [33]. We now introduce a slight inconsistency in the present formula- 
tion and assume that uniform tractions Ap (consistent with the macroscopic 
stresses AS, Ap = AS • n, n being the unit outward normal to the boundary 
of the PUC) rather than displacements are prescribed to get 



(5s : Aff(x)) = 5E : AS. 



(24) 



Substituting the discretized form of the local strain increment {Ae(x)} = 
{AE} + [B] {Ar*} together with local constitutive equation into Eq. (24) 
gives 



i/jB] T [L(x)] dfii^[B] T [L(x)] [B] dQ 




{AT}--J^{AX}dQ 
[B] T {AX}dQ 



1 

n 



(25) 

Recall that in the framework of finite element method the matrix [B] is often 
termed the geometrical matrix and the vector {Ar*} stores the increments of 
the fluctuation nodal displacements. Finally, after rewriting the above equa- 
tion as 



[K] n [K] 12 
[K] 21 [K] 22 



( {AE} ) 

1 {Ar*} J 



{AT} + {AF } 
{Af } 



(26) 



and eliminating the fluctuating displacements vector {Ar*} we arrive at the 
incremental form of the macroscopic constitutive law 



{AT} 



horn 



{AE} + {AA}, 



(27) 
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where 

|^hom 



12 i 



{AA} = -{AF°} + [K^K^Af } 



Also note that, when returning back to our original formulation with the 
prescribed overall strain only, the system of equations (26) reduces to 

[K] 22 {Ar*} = -[K] 21 {AE} + {Af }. (28) 



2.2-4 Effective properties from binary images 

A simple example of a two-scale elastic homogenization, recall Fig. [2j is pre- 
sented here for the binary images plotted in Figs. [4] and |5j Both the aggregates 
and mastic are assumed to be isotropic. The elastic modulus of the mastic 
phase being equal to 1,600 MPa [y m = 0.4) is associated with the value of 
the storage modulus at zero value of the loss modulus based on the Cole Cole 
graph in Fig. [To|(d). The elastic material properties of aggregates are consid- 
ered to be those of a spilite rock with the elastic modulus equal to 60,000 MPa 

(y 8 = 0.2). 

Table 1 

Effective properties of_ mortar matrix 



PUC (Fig. 


5 


) 


(a) 


(b) 


(c) 


(d) 


FEM (d) 


MT (d) 


E hom [ M p a j 
Z> om [-] 


2443 
0.38 


2670 
0.38 


2897 
0.38 


3150 
0.38 


3145 
0.38 


3462 
38 



Since dealing directly with the binary representation of real microstructures 
the FFT method was adopted in this study. The homogenized material prop- 
erties on the mortar scale are listed in Table [TJ Note that due to an hexagonal 
arrangement of stones within the PUC the transverse material isotropy was 
generated. In the second homogenization the effective properties of the mortar 
are then used in place of the original mastic matrix when treating individual 
RVEs in Fig. |4j The final results, promoting the proposed two-step upscal- 
ing homogenization scheme, appear as dash-pattern columns in Fig. [7j Also 
notice the checkerboard columns, derived for the same microstructures but us- 
ing mastic material properties for the binder, clearly indicating a considerable 
impact of eliminated stones on the predicted effective properties. 



2. 3 Statistically equivalent periodic unit cell 



Although considerable savings in computational time can be achieved with 
coarse RVEs, the analysis employing original microstructures still presents a 
significant challenge particularly in view of a large scale nonlinear analysis of 
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Fig. 7. Effective elastic modulus derived from images in Fig. [4] 

multi-layered rodes carried out, possibly, in the framework of multiscale FE 2 
analysis [MBS]- Even at moderate temperatures the stiffness of the binder 
might reduce by several orders of magnitude resulting in a drastic mismatch 
between the elastic moduli of binder and aggregates and consequently in a poor 
convergence performance of the FFT method. Although some improvements 
to the original formulation exist, they will not be considered in the present 
study. Instead a rather different route is built taking account of the so called 
statistically equivalent periodic unit cell [T7] . 

In particular, suppose that the original microstructure can be replaced by a 
certain artificial periodic unit cell that, from the microstructure point of view, 
statistically resembles the real material system in terms of, e.g. the two-point 
probability function, see e.g. [36] for in depth discussion on various statistical 
descriptors. Such a unit cell can be defined by the following parameters: num- 
ber of aggregates having elliptical shape, size, position, orientation and aspect 
ratio of the axes of individual ellipses. The size of stones is derived based on 
the cumulative distribution function in Fig. |3j For example, if 10 stones are 
selected for a PUC then the smallest stone corresponds to an average size of 
10% of the smallest stones determined from the cumulative distribution func- 
tion. The next stone then reflects the size of the subsequent 10% stones, etc. 
Examples of such unit cells are depicted in Fig. [8j 




(a) (b) (c) (d) 

Fig. 8. Examples of SEPUCs corresponding to a binary image in Fig. |4^d): (a) 
SEPUC 6, (b) SEPUC 37, (c) SEPUC 43, (c) SEPUC 48 

These unit cells were derived by matching the two-point probability function 
of the original microstructure, Fig. ^d), and the SEPUC. The underlying op- 



12 



timization problem was solved with the help of the evolutionary algorithm 
GRADE [37.38J . See also [T7HTB] for other applications of genetic algorithm 
based solution strategies. It is worth mentioning that no interpenetration con- 
strain was introduced as it was naturally enforced through the consistency of 
volume fractions of aggregates in the SEPUC and targeted microstructure in 
Fig- §d)- 




/ 81 1 /H l )l) 



AVE SAGE Eh,., 



Fig. 9. Effective elastic modulus derived from SEPUCs in Fig. [7] 

As typical of genetic algorithms based optimization procedures, each run re- 
sults in a unique SEPUC with a slightly different arrangement of stones, see 
Fig. |8j It therefore remains to confirm that all cells provide the "same" macro- 
scopic response. Note that for the sake of efficiency the target microstructure 
in Fig. [3|d) with mastic binder being replaced by mortar phase (the fourth col- 
umn m Table Q was used when generating artificial periodic microstructures. 
The corresponding effective properties are represented by the first column in 
Fig. [9] The remaining columns refer to the effective elastic modulus for the 
selected SEPUCs (100 such cells were generated). Although slightly differ- 
ent in their geometrical details they all provide nearly the same macroscopic 
response almost identical to the original microstructure. 

Table 2 

Effective properties of asphalt 



SEPUC (Fig. 


8 


) 


(a) 


(b) 


(c) 


(d) 


MT (d) 


Fig. 


3 


'A) 


Fig. 


1 


c) 


E hom [ M p a ] 
l/ hom [-] 


8446 
0.33 


7811 
0.35 


7934 
0.33 


7834 
0.33 


7426 
0.33 


8322 
0.32 


8998 
0.31 



The Mori-Tanaka prediction found from the same upscaling procedure as- 
suming aggregates to be represented by statistically uniform distribution of 
circular cylinders is also presented and compared in Table [2] with the selected 
estimates provided by SEPUCs (first four columns), the reduced microstruc- 
ture in Fig. [3^d) (the last but one column) and the original microstructure in 
Fig. [ljc) (last column). Although still within an acceptable range the Mori- 
Tanaka method slightly underestimates the macroscopic elastic response. This 
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difference would be even more pronounced if using mortar properties given by 
FFT simulations. 

In most practical applications, however, the mastic asphalt is typically loaded 
beyond the elastic regime. The issue, whether the geometrical invariance of 
SEPUCs outlasts even for a non-linear response should be examined. This will 
be the principal objective of the next section. 



3 Generalized Leonov model for Mastic Asphalt mixture from two- 
step homogenization 



It has been experimentally observed [39] that both asphalt binders and asphalt 
mixtures exhibit a considerable dependence on the rate of loading and tem- 
perature. This phenomenon can be well described by the Generalized Leonov 
(GL) model. Although some experimental results point out a pressure depen- 
dent behavior of these systems, the present study is limited to a Mises-type 
material with negligible plastic volumetric deformation. The essential prop- 
erties of the GL model are presented in the next section. Section 3^ then 
summarizes the experimental program to calibrate the GL model on mastic 
scale. Upscaling towards macroscopic viscous properties of the GL model is 
then addressed in Sections 13.31 - I3.4L 



3. 1 Basic properties and implementation of GL model 



Combing the Eyring flow model for the plastic component of the shear strain 
rate 

de„ 



dt 



— - sinh — . 

2A t ' 



(29) 



with the elastic shear strain rate de e / dt yields the one-dimensional Leonov 
constitutive model 140 1 



de 
dt 



de* 



+ 



de r 



de* 



+ 



r 



dt ' dt dt ' r](de p /dt) 

where the shear-dependent viscosity rj is provided by 

Vqt 



r](de p / dt) 



t sinh (t/t ) 



(30) 



(31) 



In Eq. (29), A and To are material parameters; a a that appears in Eq. (31) 
is the stress shift function with respect to the zero shear viscosity 770 (viscos- 
ity corresponding to a viscoelastic response). Clearly, the phenomenological 
representation of Eq. (30) is the Maxwell model with the variable viscosity 
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rj. Experimental evidence proving analogy between stress and temperature al- 
lows us to extend Eq. (31 ) by adding a temperature (T) dependent shift factor 
a T to get 

T)( de p / dt) = r] a a (r)a T (T). (32) 



To describe multi-dimensional behavior of the material, the generalized com- 
pressible Leonov model, equivalent to the generalized Maxwell chain model, 
can be used |41j . The viscosity term corresponding to the /x-th unit receives 
the form 

V/t = Vo,^a(r cq )a T (T), (33) 
where the equivalent shear stress r eq is provided by 




and s is the deviatoric stress tensor. Admitting only small strains and isotropic 
material, a set of constitutive equations defining the generalized compressible 
Leonov model can be written as 




^ _ V" nr ( — _ de P^ 
dt ~ h I dt dt 



a m = Ke v , (35) 

M 

s = s ^ ( 36 ) 
de 

s^ = 2 Vll ^ = 2 Vo ^a a (r cq )a T (T)^, (37) 

where a m is the mean stress, e v is the volumetric strain, K is the bulk modulus 
and is the shear modulus of the /x-th unit. 

While a fully implicit scheme has been implemented in [21] to integrate the 



system of Eqs. (35) - (37), a fully explicit Euler forward (EF) integration 
scheme is adopted here for simplicity. Thus to avoid numerical instabilities a 
relatively short time step must be prescribed. This is merely attributed to the 
fact that viscous parameters vary nonlinearly with stress, but in the context 
of EF scheme are assumed constant over the time step. 

Provided that the total strain rate is constant during integration a new state 
of stress in the matrix phase at the end of the current time step assumes the 
form 



<T n {U) = o-m(ti-i) + KAe v , (38) 
{s(ti)} = {s(ti_0} + 2G(* i _ 1 ) [Q] {Ae} + {AAft-J}, (39) 

where tj is the current time at the end of the i-th time increment; a m is the 
elastic mean stress, {s} stores the deviatoric part of the stress vector {a} and 
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{Ae} is the deviatoric part of the total strain increment. Assuming the shear 
relaxation modulus G(t) can be represented by a Dirichlet series expansion in 
the form 

G(t) = £G,exp(--^V (40) 



n=i 



and with reference to the EF integration scheme the time dependent variables 
at time instant 1V-1 receive the form 



q_S^q ^ Q cr(^~i)o-r(^~i / 
M 



At 



{AA} = -£ 

M=l L 



1 — exp 



1 — exp 
At 



At 



^a<r(*i-l)flr(*i-l) , 



At a (T (tj_i)a T (tj_ 1 ) 

{s M (^-x)}, 



, (41) 
(42) 



where represents the elastic shear modulus in the /z-th unit of the Maxwell 
chain model, 6^ is the associated relaxation time, {s M (tj_i)}, \i = 1,2, .. . , M, 
is the deviatoric stress vector in individual units evaluated at the beginning 
of the new time increment At = t{ — tj_i, and M is the assumed number of 
Maxwell units in the chain model; a CT (£j_i) and ar(ti-i) are t ne stress and 
temperature shift factors, respectively. The relaxation shear moduli G^ are 
found from a Dirichlet series expansion of creep compliance function presented, 



e.g. in Fig. 10 c) by invoking the correspondence principle. An example of a 



direct application of this principle is given in Appendix. 



According to Eq. (31) the stress shift factor a a reads 



a a (t 



i-ly 



r M £ *-l) / ginh T eg(^i-1, 



TO T 

where the equivalent stress r e(? (t i _ 1 ) follows from 



(43) 



T~eq (t 



'1 



i-lj 



(44) 



and 



[Q]=diag 



1 1 1 

1,1,1,-,-,- 
' ' '222 



(45) 



The temperature shift factor is fitted to Williams-Landel-Ferry (WLF) equa- 
tion written as 



ax = exp 



C 2 + {T-T )r 



(46) 
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where C±, C% are model parameters, To is the reference temperature and T is 
the actual temperature. 



3.2 Rate and temperature dependent behavior of mastic from laboratory ex- 
periments 



Referring to (42] we anticipated to derive the model parameter r , needed in 
turn for the determination of the shift factor a a , from a set of experiments 
when loading individual specimens in shear at constant strain rate. The avail- 
able experimental devices, however, failed to provide the required viscous flow 
at zero elastic strain increment [UJ. Additional difficulties in controlling the 
applied normal pressure when twisting the specimens immediately promoted 
the following alternative approach: 

• Construction of the compliance master curve from the measurements of 
dynamic moduli at a selected range of frequencies and temperatures to first 



get the parameters of Eq. (46) approximating the temperature shift factor 
ax- 



• Determination of the model parameter tq in Eq. (43) by horizontally shift- 
ing the experimentally derived compliance functions obtained from creep 
measurements for a given temperature at various stress levels. 

The experimental program for the determination of frequency characteristics 
of the mastic phase was carried out jointly at the Central laboratory of Eu- 
rovia Services, Ltd. and the department of Road Structures at the Faculty of 
Civil Engineering, Czech Technical University in Prague. Standard dynamic 
shear rheometer (DSR) tests were conducted under cyclic loading undergoing 
a temperature and frequency sweep from to 80°C and 0.1 to 100 Hz, respec- 
tively. For low temperature tests (from to 20°C) a mastic film was placed 
between two plates (a lower fixed and an upper oscillating) of diameter d = 8 
mm and pressed to maintain a constant height h = 2 mm during oscillations. 
For high temperatures (from 30° to 80°C) the d/h ratio equal to 25/1 mm 
was used. The examined mastic consisted of 37% of limestone filler and 62,5% 
of polybitumen AMe 65 asphalt with bulk density p m =1025 kgm -3 . The filler 
aggregates were in the following grading range: 2 mm sieve - 100%, 0.125 mm 
sieve - 81.1%, 0.063 mm sieve - 70.3%, with p s =2730 kgm" 3 . 

The storage {G ) and loss (G ') dynamic moduli were found by loading a mastic 
specimen in shear under strain controlled regime 

7 = 7sin(ta;), resulting in r = Tsm(tuj + 5), (47) 

where 7 and u are the prescribed shear strain and angular frequency, re- 
spectively, and 5 represents the phase shift between stress and strain. The 
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corresponding moduli then follow from 



G 



1 



: cos(<5), 



G 



T 



7 



sin(<5). 



(4* 



The frequency sweep was selected such as to keep the resulting stresses within 
the viscoelastic limit. 
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o 

E 
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Storage modulus G' 
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200 300 400 
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(c) (d) 
Fig. 10. (a) Dynamic moduli, (b) Shift factor, (c) Master curve, (d) Cole Cole graph 

The master curves of dynamic shear moduli were constructed by superimpos- 
ing the dynamic data measured at different temperatures for a given range of 
frequencies. Accepting the temperature superposition principle the resulting 



master curves seen in Fig. 10 a) were simply derived by horizontally shift- 
ing individual curves to comply with the data one would obtain for the 40°C 
reference temperature at much longer (data for temperatures above 40°C) 
or shorter (data for temperatures below 40°C) times if performing standard 
creep measurements. The dynamic moduli master curves were subsequently 
converted into a time-dependent compliance function at 40°C employing the 
Fourier transform. The result appears in Fig. [To](c) . For the available tem- 
perature range the resulting master curve can reliably cover the viscoelastic 
creep response over the time span up to 1000 s. To enhance numerical sta- 
bility this curve was artificially prolonged by two additional decades, see the 
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dashed line in Fig. |10[c). The associated temperature dependent variation of 



shift factor a? fitted to Eq. (46) is plotted in Fig. 10 b). The corresponding 
parameters are stored in Table |3j Finally, The experimentally-obtained values 
of G and G are plotted in the so-called Cole-Cole diagram in Fig. 10 d) used 



in the present study to estimate the elastic modulus of mastic phase, recall 
Section 2.2.4 Note that the master curves in Figs. 10 a), (c) were provided 



directly by a built in software of DSR devices. 

The second class of experiments delivered the creep data at constant temper- 
ature T = 30°C but at different magnitudes of the applied load. These data 
together with the previous results enabled us to obtain the necessary parame- 
ter r . To that end, the original master curve in Fig. [l~0|c) was first converted 
to comply with the one for 30°C. A horizontal shift of four available creep 
compliance functions finally served to estimate the parameter To, see Table |3j 
Additional creep experiments performed at temperature 50°C were used to 
validate not only the estimated value of r but also the temperature-stress 
superposition principle expressed, with reference to so called reduced time, as 



dt 



o a T [T(t)} x a a [T eq (t)}' 



(49) 



The results shown in Fig. 11 confirm applicability of the GL model to well 
represent the temperature and rate dependent behavior of mastic material. 

Table 3 



Model parameters of Eqs. (31) and (46) 



Scale 


To 


Ci 


c 2 


r [Pa] 


Mastic 


40°C 


38.56 


195.42 


1150 


Mortar 


40°C 


38.56 


195.42 


2600 


MAm 


40°C 


38.56 


195.42 


2951 



Further support for the use of GL model can be gained from the results pre- 
sented in Fig. 
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Therein, the numerical predictions of static behavior at 
different strain rates and temperatures are compared with experimentally ob- 
tained data. Note that this set of experiments was not considered to calibrate 
the material model. Realizing a difficulty of performing reproducible tests with 
such a material suggests almost a "remarkable" match. 



3.3 Response of mortar from virtual experiments 



Following on the heels of experimental work discussed in the previous section, 
a virtual set of experiments is proposed to arrive at a homogenized master 
curve and associated temperature and stress dependent shift factors on the 
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(a) (b) 

Fig. 12. Tensile loading at constant strain rate - experimental measurements vs. 
numerical predictions: (a) T = 30°C, (b) T = 50°C 



scale of mortar. To be consistent with Section 2.3 the PUC in Fig. 5[d) is 



chosen to represent an RVE for numerical analysis. The concept of first or- 



der homogenization of periodic fields outlined in Section |2.2.3| is given the 
preference to deliver the homogenized creep response of the mortar phase at 
different temperature and stress levels. 



First, a uniformly distributed range of temperatures from 0°C to 100°C was 
considered to provide for a temperature dependent viscoelastic behavior of 
mortar loaded in shear by the remote stress Ti yx = lkPa. Individual curves 
plotted in Fig. [T3~|a) were then horizontally shifted to give the homogenized 
creep compliance master curve seen in Fig. [l3^b) and consequently estimates 
of the parameters of WLF equation (46). The parameters Ci, C2 for the refer- 
ence temperature T = 40°C are stored in Table. [3] surprisingly suggesting no 
differences between individual scales. The stiffening, observed for high tem- 
peratures, is attributed to a volumetric locking owing to a very low shear 
modulus approaching to zero. This in turn yields the Poisson ratio close to 0.5 
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Fig. 13. (a) Creep data at different temperatures for reference stress T, xy 
(b) Master curve for reference temperature T = 40°C 



lkPa, 



in a finite zone of the binder phase that is progressed over the entire unit cell 
as seen in Fig. [l3^b) . Note that such a microstructure dependent response can 
hardly be represented by the Mori-Tanaka method, see appendix for further 
illustration. 



The second set of creep experiments was conducted at two different tempera- 
tures and two different levels of the remote stress H xy . The two representative 



results, [40°C, lOkPa] and [0°C, 20kPa], are plotted as solid lines in Fig. [13(b). 
The indicated horizontal shifts identified with the corresponding star-lines 
then supplement the necessary data for the derivation of stress dependent 
shift factor a a fitted to Eq. (43). The searched parameter r is available in 
Table. K3) Note that prior to shifting the curves the result derived for 0°C was 



thermally adjusted to be consistent with the 40 C master curve. 



3.4 Response of M Am from virtual experiments 



Derivation of the macroscopic creep compliance master curve for a Mastic 
Asphalt mixture follows the general scheme sketched in the previous section. 
A statistically equivalent periodic unit cell introduced in Section 2.3 is con- 
sidered for numerical simulations. The principal assumption of restricting the 
material model on every scale to one particular format is brought to light here 
by suggesting the binder be well represented by the homogenized GL model 
predicted on the mortar scale. 

However, for problems involving significant nonlinearities, the fundamental 
statement borrowed from elasticity solutions that any statistically equivalent 
periodic unit cell is equally suitable should be scrutinized first. 

To address this issue the statistically equivalent unit cells in Fig. [8] were sub- 



21 




(a) (b) 

Fig. 14. (a) Macroscopic response for various SEPUCs T = 20°C, E xy = 10~ 4 , (b) 
Master curves on individual scales for reference temperature T = 40° C 



jected to a remote shear strain rate E xy = 10 . The resulting homogenized 



stress-strain curves for temperature T = 20°C are shown in Fig. |l4|(a). Al- 
though no "perfect match" is observed, the difference in estimated load bear- 
ing capacities is in the same range as the predicted elastic constants, see e.g. 
Table |2j not exceeding 10%. On the contrary, the distribution of local fields 



varies considerably as seen in Fig. 15 While a highly localized distribution 



of shear strain 7™ in the mortar phase, identified with the lowest bearing 
capacity in Fig. |l4|(a), is evident in Fig. [i~5^a), the variation of this quantity 
in Fig. 15 c) shows a rather distributed character consequently resulting in a 
slightly stiffer response on the macroscale. 



Although certainly more accurate, the detailed finite element simulations are 
in general computationally very expensive and often call for less demanding al- 
ternatives such as the Mori-Tanaka method. Unfortunately, the corresponding 
results also plotted in Fig. [l~4"|(a) clearly expose essential limitations of two- 
point averaging schemes, unable to capture localization phenomena observed 
in composites with a highly nonlinear response of the binder phase [21J. On 
that ground, a considerably more refined variant of the TEA method would 
be necessary |13J3I]. In this paper, however, the attention is limited two a 
two-phase composite only as discussed in Section [2.2. 1| The presented results 
in Fig. [l4|a) correspond to a classical formulation with the localization and 
transformation tensors calculated only once being functions of elastic prop- 
erties of individual constituents (MT-original), and to the formulation where 
these tensors are updated after each time step taking into account increasing 
compliance of the mortar phase with time (MT-improved). Note that even 
the latter case produces much stiffer response, although at a fraction of time, 
when compared to the FE results. 



Regarding the "similarity" of macroscopic response from various SEPUCs, the 
SEPUC No. 37 was selected to provide data needed in the calibration of the 
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(b) 



(c) 



(d) 



Fig. 15. Distribution of local shear strain: (a) SEPUC 6, (b) SEPUC 37, (c) SEPUC 
43, (d) SEPUC 48 



macroscopic GL model. Virtual numerical tests identical to those in Section [373 
were again performed to give first the homogenized master curve displayed 



in Fig. 14 b) and consequently the model parameters of Eqs. (43) and (46) 



available in Table |3j Note again the same temperature shift as obtained already 
for lower scales. This is evident in Fig. |14~[b) showing only vertical shift of 
individual master curves attributed to an increasing stiffness when moving up 
individual scales. 



Finally, fitting the inverse of the homogenized master curve to a particular 
chain of Maxwell units then allows for estimating the macroscopic response 
of MAm to a given load limiting the material symmetry to a macroscopic 
isotropy. The response of the homogenized asphalt mixture to the applied 
remote shear strain labeled as "Leonov macro" appears in Fig. |l4[a). Since 
derived from the solution of a unit cell problem, a relatively good agreement 
with these solutions has been expected. The difference is merely attributed 
to the specific approximation of the homogenized master curve via Dirichlet 
series. It is fair to mention a considerable sensitivity of the predictions to 
a particular choice of pairs of parameters (shear modulus and retardation 
time) of the Maxwell units. Although promising, to fully accept the proposed 
uncoupled multiscale computational strategy will require further testing, both 
numerical and laboratory. This topic is still widely opened. 



4 Conclusions 



Although research interests on flexible pavements have been quite intense in 
the past two decades, the field is still very much in development and will cer- 
tainly witness considerable activity in the coming decade particular in connec- 
tion to hierarchical modeling and micromechanics. Within this framework, the 
present contribution provides theoretical tools for the formulation of macro- 
scopic constitutive law reflecting the confluence of threads coming from exper- 
imental work, image analysis, statistical mechanics and traditional disciplines 
of micromechanics and the first order computational homogenization. Here, 
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the totally uncoupled multiscale modeling approach is favored to enable an 
inexpensive analysis of real world large scale structures. Since much of the 
considered is primarily computational it would be audacious to brought this 
approach directly to points of application without additional experimental 
validation. Large scale experiments of rut depth measurements due to moving 
wheel load are currently under way. The effect of microstructure anisotropy, 
influence of the first stress invariant or three-dimensional character of asphalt 
mixtures are just a few issues which need to be addressed, yielding flexible 
pavements a fertile field of future research. 
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A The Mori-Tanaka predictions based on correspondence principle 



The correspondence principle as a tool allowing for interconversion of linear 



viscoelastic response functions has already been mention in Section 3.1 Al- 
though not used in the present contribution, this principle is reviewed here 
to fill in the list of possible routes for deriving the homogenized viscoelastic 
response of composites using simple averaging schemes. As an example, the 
Mori-Tanaka method is adopted again to construct the homogenized master 
curve of a two-phase composite medium composed of aligned elastic circular 
fibers bonded to a viscoelastic mastic phase. This approach thus allows direct 
comparison with the finite element predictions obtained on the mortar scale. 

Consider a homogenized relaxation loading path written as 



M / 

G hom (t) = ^Gj om exp - 



flhom^hom^hom I ' 

> "cr a T ) 



(A.l) 
(A.2) 



where M is the number of Maxwell units representing the homogenized re- 
,iiom — i [ n the viscoelastic regime, a)j? m = 1 for the reference tem- 
40°C (assumptions complying with the corresponding master 
13 b)) and 9 h ° m are appropriate relaxation times. Then, 



sponse, „ a 
perature T 
curve shown in Fig. 
according to the correspondence principle the Laplace-Carson transform of 



Eq. ( A.l ) gives 
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K y {v) = Gl om {j))E* xy {jp), 
M G horn 



G* hom (p) = E 



' Ohom 

ft 



+ P) 



(A.3) 
(A.4) 



Following [15] provides Eq. (12) in the Laplace space in the form 



m 



MT 



m *s m m(.krn + ^ m rn) + Kn m *rn{ C s m *s + c m m m 

/c^m^ + (k* m + 2m* m )(c s m* m + c m m*) 



(A.5) 



Since the volumetric response is assumed to remain elastic, recall Eq. (35), we 
get 



Kip) 



m: 



p(l - 2i/„ 



M=l ( h n) 



m s 

P 

V 



(A.6) 
(A.7) 

(A.8) 



where N is the number of Maxwell units approximating the behavior of the 
mastic phase and a™ = a™ = 1 is again adopted to guarantee compatibility 



with the results presented in Section [372] Once the Mori-Tanaka estimate Gmt 
of the homogenized relaxation function Ghom in the Laplace space are known, 



the desired coefficients G^ om , which enter Eq. (A. 2), are found by solving the 
following minimization problem 



K 



e g = J2 l G *Mj(Pk) - G* hom (pk)} = min, 



(A.9) 



thereby avoiding a cumbersome and rather unstable numerical transformation 



of Eq. (A.5) back to time dependent solution (A. 2). Solution strategies based 



on soft computing [37|35 ,38j are typically employed to solve Eq. (A.9). 



Given the relationship between relaxation and compliance functions 



pJ*ip) 



pG* (p) 



(A.10) 
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allows us to write a dual representation of Eq. (A. 5) in the form 



P Jmt(p) 



m 



MT 



+ k m + 



P 2J m) \P 2 Jn 



+ c m m a 



P 2 J* 



km + 



+ 



k* 



c s m c + 



(A.ll) 



where 



N 



J*miP) = 



jm 



(A.12) 



The homogenized coefficients Jfi hom then follow from the minimization of a 
dual error function 



K 



mm. 



(A.13) 



Pk 



or directly from Eq. (A. 10) providing the set of coefficients Gfi hom has been 



already found from the primary formulation. 



The resulting MT prediction of Jmt for the reference temperature T = 40°C 
is compared in Fig. A.l[ b) with the PUC solution on the mortar scale, recall 
Fig. 13 b). Inability of the MT method to capture the stiffening effect already 



mentioned in Section 3^ is evident in both the Laplace and time space. It is 
interesting to point out that such a discrepancy appears already for the elastic 
predictions if setting m m — > and v m — > 0.5, in which case the binder behaves 
as an incompressible fluid. 
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